Cumulative distribution functions associated with bubble-nucleation processes in 

cavitation 
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Bubble-nucleation processes of a Lennard- Jones liquid are studied by molecular dynamics simu- 
lations. Waiting time, which is the lifetime of a superheated liquid, is determined for several system 
sizes, and the apparent finite-size effect of the nucleation rate is observed. From the cumulative 
distribution function of the nucleation events, the bubble-nucleation process is found to be not a 
simple Poisson process but a Poisson process with an additional relaxation time. The parameters 
of the exponential distribution associated with the process are determined by taking the relaxation 
time into account, and the apparent finite-size effect is removed. These results imply that the use of 
the arithmetic mean of the waiting time until a bubble grows to the critical size leads to an incorrect 
estimation of the nucleation rate. 
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I. INTRODUCTION 



When a system exhibits a first-order transition at some 
set of parameter values, nucleation phenomena are ob- 
served in a transient process. A familiar example of such 
nucleation is observed in the liquid-gas phase transition. 
The theory describing such nucleation was initiated by 
Gibbs [l] , which was followed by quantitative arguments 
for the steady-state nucleation rate @,y]. This nucleation 
theory was refined by Zeldovich [jand is now called the 
classical nucleation theory (CNT) [51]. CNT was first con- 
structed for droplet nuclei in a supersaturated vapor and 
then applied to bubble nuclei in a superheated liquid [g] . 
In CNT, the nucleation is treated as a stochastic process, 
and CNT predicts a nucleation rate that corresponds to 
the emerging frequency of critical embryos. While nu- 
cleation rates during the homogeneous condensation of a 
supersaturated vapor are, on the whole, well predicted by 
CNT with some modifications [7], it is well known that 
the nucleation rates of bubbles in a superheated liquid 
predicted by CNT are markedly different from the val- 
ues obtained from experimentally [8i] or numerically [9|. 
This discrepancy originates from the fact that the nu- 
cleation phenomena of bubbles are considerably different 
from those of droplets for the following reasons, (i) The 
free energy required to form a bubble is not a unique 
function of its volume or the number of particles since a 
gas is compressible, (ii) Work carried out by bubbles on 
the ambient liquid cannot be ignored, (iii) Interbubble 
interactions via ambient liquid cannot be ignored. One of 
the important differences between the droplet and bub- 
ble nucleation is the density of the ambient phase. For 



the case of the droplet-nucleation, the ambient phase is 
gas phase which density is usually negligible. However, 
the bubbles can interact each other via ambient liquid. 
For example, one bubble growth increases the pressure of 
the surrounding liquid which may suppress the nuclei of 
other embryos. Corresponding to (i), the compressibil- 
ity of bubbles can be measured by molecular dynamics 
(MD) simulations [!£], and the free energy surface of the 
embryos of bubbles has recently been studied in terms of 
density functional theory (DFT) [Tl|. The work carried 
out by bubbles can also be estimated directly by MDja] 
and indirectly by DFT from the analysis of cavity [l2l [. 
Despite such studies, less attention has been paid to the 
finite-size effect on bubble-nucleation, while such effect 
was investigated for droplet-nucleation [13| . If the effect 
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from (iii) is exhibited, this would cause strong finite-size 
effect. Therefore, investigation of the size dependence of 
nucleation rates is necessary. In the present article, we 
study bubble-nucleation processes of a Lennard- Jones liq- 
uid by MD to clarify the finite-size effect of bubble nuclei. 



II. NUCLEATION RATE 

Consider a superheated liquid in a metastable state. 
While the gas phase is more stable than the liquid phase, 
time is required for the uniform phase to change to the 
gas-liquid coexistent phase since there is an energy bar- 
rier to be overcome in the formation of large bubbles. The 
lifetime of such a superheated liquid, which we call the 
waiting time t^ in the following, is a stochastic variable. 
The superheated liquid is characterized by a nucleation 
rate J, which is the number of embryos growing beyond 
the critical size in the unit time and volume. If the sys- 
tem is in the steady state, J can be expressed in terms 
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FIG. 2: Phase diagram showing the vicinity of the phase 
boundary between the pure hquid and gas-Uquid coexistent 
regions. The solid hne denotes the binodal hne between the 
liquid and gas-liquid coexistent phases, and the dashed line 
denotes the spinodal line where spinodal decomposition takes 
place. The filled and open squares on the binodal and spin- 
odal lines denote the observed point which are used to draw 
the lines. The system is first thermalized in the pure-liquid 
phase, then is suddenly expanded. Then the system crosses 
the binodal line and becomes a metastable state (superheated 
liquid) in the vicinity of the spinodal line. The solid and open 
circles denote the states of the system before and after expan- 
sion, respectively. (Inset) Overall view of the gas-liquid phase 
diagram. The symbols 'G', 'L', and 'GL' denote the pure gas, 
pure liquid, and gas-liquid coexistent regions, respectively. 



FIG. I: Time evolutions of (a) temperature and (b) pressure. 
Only the data for L — 64, 128, and 256 are shown for visibility. 
After thermalization at T = 0.9, density suddenly decreases 
from 0.7 to 0.659 at t = 50. Then the system relaxes to a 
metastable state with T — 0.868(2). The time evolutions of 
different system sizes are almost identical. 



of the expectation value of the waiting time (iw) as 

1 



J = 



L^ (iw) ' 



(1) 



with the linear size of the system L. It is widely assumed 
that a bubble-nucleation process is a Poisson process, 
then the cumulative distribution function (CDF) associ- 
ated with nucleation events has the exponential distribu- 
tion 



F{t) = P(tw < t) == 1 - exp(-i/r). 



(2) 



The exponential distribution is specified only by the pa- 
rameter T, which denotes the time scale of this stochastic 
process. From Eq. ([2]), {t^j) equals r. 

On the other hand, CNT predicts the nucleation rate 
from the work carried out to form the critical bubble as 



J = DZc{n*,v*), 



(3) 



where D is the kinetic pref actor, Z is the Zeldovich fac- 
tor, which describes the nonequilibrium effect, and c{n, v) 
is the number density of bubbles with volume v contain- 
ing n particles [6|. In a superheated liquid, the reversible 
work W(n,v) carried out to form a bubble with volume 
V containing n particles has a saddle point at {n*,v*), 
which corresponds to the critical bubble. The nucleation 
rate J given in Eq. (j3|) is determined only by the temper- 
ature and density of the superheated liquid. Equations 
(H)) and ([3]) lead to the finding that (iw) should be in- 
versely proportional to the volume of the system if the 
system does not have any finite-size effect. 



III. METHOD 

To estimate the waiting time in bubble nucleation, we 
perform MD simulations. We use the truncated Lennard- 
Jones potential of the form 



V{r) = 4e 



(7) 
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cr\" /r 

- +C2 - 
r / \a 



Co 



(4) 



with the well depth e and atomic diameter a [M]. The 
coefficients C2 and cq are determined so that V(rc) = 
V'i'i^c) = with the cut-off length re , i.e., the values 




FIG. 3: (Color online) Time evolution of bubble in the system with L = 128. Left to right: snapshots at f = 100, 150, and 
550. We divide the system into small subcells and define them to be in the gas state when their densities are less than 0.2. 
Only the gas state subcells are shown. An embryo appears after a waiting time then grows slowly into a large bubble. 
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FIG. 4: System-size dependence of the expectation value of 
the waiting time (iw). Decimal logarithms are taken for both 
axes. The error bars are smaller than the symbols. The solid 
line, which denotes L~^ , is a guide for the eyes. 



of potential and force become continuously zero at the 
truncation point. In the following, we use the physical 
quantities reduced by cr, e, and k^, i.e., the length scale 
is measured using the unit of ct, and so forth. We set 
the cutoff length as Tc — 3.0. The system is a cube with 
linear size L and is periodic in all directions. The num- 
ber of particles is chosen so that the initial density of the 
system p = N/L^ is 0.7. We first maintain the system in 
the pure-liquid phase using a thermostat, then we expand 
the system. The expansion is performed by changing the 
radius of the particles from a to a' — aa, where a is a 
rescaling factor. This procedure is equivalent to the the 
uniform and adiabatic expansion, which is q.j — >■ Qj/a 
and L —> L/a where q, is the position of the particle i. 
Note that all physical quantities should be rescaled after 



expansion, since we measure them in the unit of the ra- 
dius (7, for example, T' = T/a^, p' = po? ^ and so forth. 
We chose the rescaling factor a to be 0.98 for all runs, 
and therefore, the change in the density upon the expan- 
sion is from p = 0.7 to 0.659 [1^. After the expansion, 
we turn off the thermostat and continue the microcanon- 
ical simulation. The system is thermalized at the tem- 
perature T = 0.9 using the Nose-Hoover method [iq . 
The integration scheme for the isothermal time evolu- 
tion is the second-order Reversible System Propagator 
Algorithm (RESPA). J7;], and the leapfrog algorithm 
is used for the microcanonical simulation with the time 
step At = 0.005. The typical time evolutions of tem- 
perature and pressure are shown in Fig. [T] One can 
see that both temperature and pressure suddenly drop 
when the systems are expanded, and they relax to values 
for metastable states, which are superheated liquids with 
negative pressure. The time evolutions of temperature 
and pressure are less affected by the size of the system. 
All systems become superheated liquids at temperature 
T = 0.868(2). We obtain the phase diagram of this sys- 
tem from the preliminary simulations. We also determine 
the spinodal line between the liquid and liquid-gas coex- 



istent phases by the method described in Ref. [18[. The 
obtained phase diagram is shown in Fig. [2l As shown in 
the figure, the system with p = 0.659 and T = 0.868 is 
in the liquid-gas coexistent region. The densities at the 
binodal and spinodal points for T — 0.868 are 0.695(1) 
and 0.640(1), respectively. In order to identify bubbles, 
we divide the system into small subcells with length 2.0 
and observe the local density for each subcell. From the 
preliminary simulations, the densities of gas and liquid 
coexisting in this system at T = 0.9 are estimated to be 
0.04(1) and 0.67(2), respectively. Therefore, we define a 
subcell to be in the gas state when its density is less than 
0.2. We have confirmed that the results do not change for 
other values of the threshold such as 0.3. We define that 



>- 



256 : 


■192 


/ 128 














. 96 






: \ 




^,..--"^" 


^,..--'"' 


64 ^,. 


32 



50 



100 
t 



150 



200 



10 



10 



•^10 



10 



10' 



3 












2 




0\ 


"v. 






1 










"\ ■ 













\,^^ 



32 64 96 128 192 256 

L 



FIG. 5: Cumulative distribution functions of nucleation 
events. The values of Y{t) = — ln(l — F) are shown as func- 
tions of time. The labels near lines denote system sizes. 



FIG. 6: System-size dependence of the characteristic time 
scale r of the Poisson process. Decimal logarithms are taken 
for both axes. The error bars are smaller than the symbols. 
The solid line, which denotes L~^, is a guide for the eyes. 



the neighboring gas state cells are in the same cluster and 
identify the bubble using the site-percolation criterion in 
the simple cubic lattice. The time evolution of a bubble 
identified by the above method is shown in Fig. [3] 



IV. RESULTS 

We first estimate the critical size of a bubble. After 
expansion, the volume of the largest bubble fluctuates 
for a certain period of time, and then the monotonic de- 
velopment of a bubble with a different waiting time is 
observed. If a bubble exceeds some critical size, then it 
starts to grow explosively. Conversely, the sizes of bub- 
bles before the explosive growth should be smaller than 
the critical size. The maximum volume of the bubbles in 
the region of fluctuation is estimated to be about v — 180. 
We therefore define the waiting time t^ as the interval 
between the time of expanding operation (i = 50) and 
the time when the volume of the largest bubble reaches 
V — 200. We study several system sizes from L = 32 to 
256. The sizes of the studied systems and the number of 
particles are listed in Table HI We observe 256 indepen- 
dent samples of the waiting time for each system size and 
take their simple arithmetic mean. Computations are 
mainly performed on HITACHI SR16000/L1 (32 ways on 
1 node). For the largest systems with 11744051 particles, 
40850 steps including those for thermalization are calcu- 
lated in an average of 12274, which give the calculation 
speed of 39.1 million updates per second. The system- 
size dependence of the waiting time is shown in Fig. 21 It 
is apparent from the figure that the waiting time is not 
proportional to L~^ , and therefore the nucleation rate 
exhibits a strong finite-size effect. 

In order to clarify the reason for the finite-size effect, 
we observe the probability distribution of the nucleation 



events. First, we determine whether the nucleation pro- 
cess is a Poisson process. Equation ([2]) leads to 



</t = -ln[l-F(i)] 



(5) 



which means that the function Y{t) = — In [1 — F{t)] be- 
comes linear and its slope gives the parameter of the dis- 
tribution provided the nucleation is a Poisson process. 
The values oiY(t) are shown in Fig.[5j Whereas the lines 
are straight for small systems, those of larger systems 
bend at low values of t. Additionally, the x-intercepts 
are not located at the origin, which suggests that the 
distribution of the waiting time has the form 



F{t) = 



■exp[-(t-fo)/T] 



t<to, 

t>to, 



(6) 



with the additional relaxation time to. Therefore, we per- 
form a fitting assuming the form given by Eq. ^. For 
systems with L > 96, we only apply a fit to the region 
where the line appears to be straight in Fig. [5] The fit- 
ting results are summarized in Table ID It is shown that 
the additional relaxation time to is almost independent of 
the system size. The additional relaxation time was also 
reported in droplet-nucleation [13, 19]. While it was con- 
sidered to be the time that the system needs to produce a 
nucleated cluster in droplet-nucleation, it can be a result 
due to the expansion, since the relaxation time does not 
exhibits the finite-size effect such as the temperature of 
the liquid under expansion as shown in Fig. [1] In order 
to estimate the relaxation time due to the expansion, we 
observed the autocorrelation function of the temperature 
after the expansion, that is, the autocorrelation function 
of the superheated liquid. We assumed that the autocor- 
relation function's form has a simple exponential form, 
C{t) r^ exp(— t/ra), with the characteristic time scale Ta, 
and we found that Ta = 2.00(2), which is too short to ex- 
plain the value of additional relaxation time to ~ 30. We 



L 32 64 96 128 192 256 

iV 22937 183500 619315 1468006 4954521 11744051 

(iw> 2.0(3) X 10^ 150(12) 68(3) 54(2) 37.9(9) 31.8(6) 

r 1854(1) 114.3(5) 32.5(3) 19.7(1) 7.63(8) 4.2(1) 

to 27(4) 36.6(2) 39.7(2) 36.3(2) 32.4(1) 29.6(1) 



TABLE I: Summary of the physical quantities for systems 
with length L and number of particles TV. The arithmetic 
mean of the waiting time is denoted by (fw). The character- 
istic time of the exponential distribution r and the additional 
relaxation time io are determined by the CDFs of iw 




also check the definition of the waiting time tw The value 
of the waiting time depends on the volume of the criti- 
cal bubble where we choose v = 200 for the nucleation 
threshold in this study. However, the value of r should 
be independent of this definition since it is a parameter 
which depends only on the state of the superheated liq- 
uid. In order to confirm the above, we observe the wait- 
ing time with different values of the threshold v — 180 
and 220 for L = 64. Then we find that the additional 
waiting time io becomes longer for larger values of the 
threshold as ig = 35.0,36.6, and 37.4 for v = 180,200, 
and 220 while the timescale parameter r is almost in- 
dependent of the value of the threshold. Therefore, we 
conclude that the additional relaxation time in bubble- 
nucleation is also a time to produce a nucleated cluster 
as in droplet-nucleation. 

The system-size dependence of the parameter of the 
exponential distribution r is shown in Fig. [51 which shows 
that the parameter is almost proportional to L~^ , which 
implies that the nucleation rate is almost independent of 
the system size. Assuming the relation J = l/(rL^), as 
suggested by Fig. [51 we estimate the nucleation rate to 
be J = 1.66(5) X 10~^, which is much larger than the 
value of J ~ 2.63 x 10"" predicted by CNT, as weU 
known [20|. The nucleation rate still exhibits the finite- 
size effect for larger systems, for example, the value of r 
when L = 256 is only five times shorter than that when 
L = 128, which should be eight times shorter without 
the finite-size effect. In order to illustrate how close the 
bubble nucleation is to a Poisson process, we plot the 
CDF as a function of the value X — 1 — exp [— (t — to)/'''] 
in Fig. [71 If the nucleation process is a stochastic process 
with the distribution given by Eq. ([S]), then the plot will 
become a straight line connecting the origin to (1,1). 
The figure shows that the CDFs of the larger systems 
deviate from the Poisson process for small values of X. 
This implies that the relaxation process caused by the 
expansion may affect bubble nuclei, but it is difficult to 
separate the time scale of bubble nuclei from that due to 
the expansion when two time scales are similar to each 
other. 



FIG. 7: Cumulative distribution function F{t) plotted as a 
function of X = 1 — exp [— (i — to)/T] for system sizes with the 
values listed in Table. [J If a nucleation process is perfectly 
described by Eq. ((5| , then this plot becomes a line connecting 
the origin to (1, 1). 



V. SUMMARY AND DISCUSSION 

In the present study, we investigated the bubble- 
nucleation process of a Lennard- Jones liquid by MD sim- 
ulations. We found that the nucleation rate, defined 
by the arithmetic mean of the waiting time, exhibits 
an apparent finite-size dependence that can be removed 
by analysis of the CDF associated with the nucleation 
events. The obtained CDF is found to be an exponential 
distribution with an extra delay time. Therefore, it is 
necessary to study a system in which the waiting time 
is sufficiently longer than the additional relaxation time 
to study bubble nuclei in a superheated liquid. The nu- 
cleation rate correctly determined by CDF has a smaller 
finite-size effect. This suggests that the discrepancy be- 
tween the prediction of CNT and the experimental data 
does not originate from the result of interbubble interac- 
tions though the pressure of the ambient liquid, but from 
the inaccurate estimation of the reversible work carried 
out to form a critical bubble. To verify this conjecture, 
the direct measurement of the reversible work is required, 
which can be achieved by observing the bubble distribu- 
tion in a superheated liquid. This is one of important 
issues. We investigated the origin of the additional re- 
laxation time to. Similar to the droplet-nucleation case, 
to is the time required to make the critical bubble from 
the candidates which are the fiuctuating embryos in the 
superheated liquid. We also studied some different ex- 
pansion ratios and found that the additional relaxation 
time hardly depends on the expansion ratio which implies 
that to depends only on the density and the temperature 
of the meta-stable liquid after the expansion. Figure [3 
shows that there are bubble nuclei in the short-time re- 
gion t < to, which is apparent for larger systems, while 



we assumed that the probabiUty of nucleation is zero in 
Eq. ^. This suggests that the inhomogeneity of density 
caused by the expansion of the system may enhance bub- 
ble nuclei, but further studies are required to clarify the 
effect of expansion on the formation of bubbles. 
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